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Abstract 
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*  For  several  years  the  authors  have  been 
involved  In  evaluating  and  developing  numerical 
simulation  schemes  for  compressible,  two 
dimensional  inviscid  or  viscous  flows  in 
turbomachinery  cascades.  Numerical  schemes 
considered  Were  all  originally  classified  as 
time-marching  schemes. and  include:  1)  implicit 
approximate  factorization  schemes- like  those  of 
Beam  and  Wanning,  2)  explicit  schemes  due  to 
MacCormackT  3)  explicit  central  difference  schemes* 
due  to  Jameson,  Schmidt  and  Turkel  and  Rizzi,  and 
4)  the  multi-grid  scheme  of  Ni,  As  we  have 
developed  these  schemes  we  have  cipme  to  believe 
that  the  accuracy  of  computational  results  is 
relatively  insensitive  to  the  numerical  algorithm 
chosen  but  highly  sensitive  to  implementation 
details  such  as  boundary  conditions,  consistent 
flux  balancing,  gri.d  resolution  and  numerical 
smoothing.  In  order  to  illustrate  our  viewpoint, 
we  present  an  examination  of  the  relationship 
between  a  flux  balancing  interpretation  of  the 
control  volume  conservation  laws  and  various  finite 
difference  formulations  and  comparisons  of  the 
performance  of  these  schemes  on  three  test 
problems:  Ni’s  bump  in  a  channel,  a  supersonic 

nozzle,  and  flow  in  a  supercritical  compressor 
cascade.  As  a  final  test  example,  a  comparison  is 
made  between  measured  and  computed  flow,  using  a 
compressible,  Navier-Stokes  solver,  for  a  high 
speed  turbine  cascade. 


Finite  Difference  Ope r ator  Definitions 

4j  uj,k  -  2<uj*1,k  '  u3-1,k> 

“j  uj.k  *  I(uj+1.k  +  uj,k> 

Uj  uj,k  *  2<uj.k  + 

*j  uj,k  *  uj+1,k  ‘  uj,k 

4j  uj,k  *  uj.K  *  uj-1.k 
4  j/2  uj.k  *  uj^,k~  uj-i.k 


Note  that  all  coordinate  positions  at  non-integer 
mesh  spacing  are  to  be  defined  by  simple  averages, 
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For  several  years  the  authors  have  been 
involved  in  evaluating  and  developing  numerical 
methods  for  compressible,  inviscid  or  viscous, 
internal  flow  problems.  The  methods  discussed  in 
this  paper  include  (1)  implicit  approximation 
factorization  schemes  due  to  Beam  and  Warming  [1), 
(2)  explicit  schemes  due  to  MacCormack  (2),  (3)  an 
explicit  central  difference  scheme  suggested  by 
Jameson,  Schmidt  and  Turkel  [3]  and  Rizzi  [4],  and 
(4)  the  multi-grid  scheme  due  to  Ni  15] .  Most  of 
these  methods  were  originally  formulated  as 
time-marching  methods  although  all  are  now  run  as 
pseudo-time-marching  methods. 


Throughout  these  evaluations  we  have  been 
interested  primari ly  in  the  accuracy  of  the  steady 
state  solutions  rather  than  the  computer  run  time 
or  number  of  iterations  to  convergence.  As  we  have 
developed  these  methods  we  have  come  to  believe 
that  the  accuracy  of  the  computational  results  is 
relatively  insensitive  to  the  numerical  scheme 
chosen,  but  it  is  highly  sensitive  to  implemen¬ 
tation  details  such  as  boundary  conditions, 
consistent  flux  balancing,  grid  resolution  and 
numerical  smoothing. 


In  order  to  illustrate  our  viewpoint  we  first 
present  an  examination  of  the  relationship  between 
a  flux  balancing  or  finite  volume  interpretation  of 
the  control  volume  conservation  laws  and  a  finite 
difference  formulation  in  strong  conservation  law 
form,  while  the  development  here  is  perhaps  not 
unique  or  entirely  new,  it  is  essential  for  under¬ 
standing  the  relationships  between  various 
numerical  schemes.  Next  we  present  a  development, 
from  a  flux  balance  viewpoint,  of  the  several 
schemes  which  emphasizes  their  basic  similarity 
rather  than  their  computational  differences.  The 
conclusion  here  is  that  all  these  schemes  should 
produce  nearly  the  same  steady  state  solutions  when 
the  same  boundary  conditions  and  smoothing 
operators  are  applied.  Results  from  three  inviscid 
computational  test  problems  are  then  presented  in 
order  to  demonstrate  the  validity  of  this 
conclusion.  Lastly,  results  from  computations 
using  a  full  Navier-Stokes  version  of  the 
Beam-Warming  algorithm  are  compared  to  experimental 
results  for  a  high  speed,  high  turning  turbine 
cascade.  These  results  illustrate  that  a  consis¬ 
tent  application  of  the  flux  balancing  analysis  can 
produce  accurate  solutions  even  in  extreme 
geometries  with  rather  modest  cost.  The  number  of 
iterations  to  full  convergence  is  about  500  for  a 
grid  containing  100*50  nodes. 
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Relationship  Between  Finite  Difference  and 
Finite  Volume  Interpretations 

Consider  the  two-dimensional  Euler  equation 
expressed  in  conservation  law  form: 


3U  3F  3G 

at  +  8x  3y  * 


pu 

DU2  +  P 
puv 

pu<E+P/p) 


pv 

puv 

pv2  +  P 
pv(E+P/p) 


and  a  coordinate  system  transformation# 

C  *  C(x,y)  n  «  n(x,y)  t  *  t  . 

The  strong  conservation  law  form  under  this 
transformation  is  retained  as#  reference  [6}: 


where  J  is  the  transformation  Jacobian 


_  (ill  (ini 

^ 3 xj  ^3yJ  ^ 3y  I  |^3xJ 


In  shorthand  notation#  equation  (3)  becomes: 

U,t  *  +  G#n  ~  0  ' 


0  *  ”  t 
J 


G«i-nx+~nv  * 
j  ,x  j  ,y 

J  “  S,x  n,y  “  ^,y  '•.x 

Formal  relationships  exist  between  the 
transformation  from  (x#y)  coordinates  to  (C,n) 
coordinates  and  the  inverse  transformation  from 
(C#n)  coordinates  to  (x,y)  coordinates  which  are: 


x  -  n  /J  # 

»y 


-  «,y/j  . 


-  \x/J  ' 


C  /J  , 


m  +  _ifr.  n 

+  S  3 k) 

+  _3(F  3n 

♦  s  -o 

[ 1 

|jJ  3C^J  3x 

J  3y  J 

3n[j  3 x 

J8yJ  (3, 

I 

and  the  inverse  Jacobian  is 

J_1  -  x.<  y.->  -  x.n  y,C  -  VJ  . 

The  Jacobians  express  the  ratio  of  cell  areas 
between  the  different  representations. 

The  transformation  metrics  (t#x  C#y  n#x 
n^y)  are  most  conveniently  obtained#  as 
suggested  by  Steger  (7) ,  from  a  finite  difference 
evaluation  of  these  derivatives  in  computational 
space  using  equations  (6).  This  process  uses  only 
the  fact  that  the  coordinate  transformation  is  a 
one-to-one  mapping.  If  we  consider  all  these 
derivatives  to  be  evaluated  by  centered  finite 
differences#  we  have,  assuming  A£  *  An  “  1# 


yjfk  - 


(Xj#k+1  ~~  xj#k-l) 
2 


(V)j,k  -  -  ‘K  *j.k  -  :{xi aa=i> 
-  -  «!  y^  - 

(v)j,k  -  6°  -  (X^'K ;  x-3',(k) 


6j  xj»k 


*x j  +  1 ,k  ~  xi-1 #k l 
2 


If  we  consider  the  computational  cells  to  be 
polygons,  the  exact  area  of  a  cell  can  be 
determined,  see  [8] ,  from 

A  N  N  1 

Area  “  L  xi  yi«1  -  £  yi  | 

where  N  is  the  number  of  edge  nodes,  and 


if  i+1  g  N 
if  i+1  >  N 


S’!**-/) 

7  *  7 (x,y) 


y*y($,’>> 


Figure  1.  Relationship  between  transformation 
from  physical  to  computational  space 
and  computational  to  physical  space. 


Figure  2.  Node  connection  diagram  for  general 
computational  cell  polygon. 
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Now  consider  a  possible  finite  difference 
representation  of  equation  (5), 


( au/j]n 

l Jj.» 


6j  rj,k  -  «k  cj,k 


on  a  standard  9  point  cell, 


j-1 ,k+1 

j  A+1 

j+ 1 , k+ 1 

j-1  ,k 

j  A 

j+1  ,k 

j-1  A-1 

j,k-1 

j+ 1 #k- 1 

Figure  3.  Finite  difference  node  numbering  scheme. 


(->  <U/J)]n  «  .  1  { (F  L*)  +  fG 

lJt  J  j,k  2  I l  J  j  j+1,k  l  J  J  j+1,> 


fr  VI 

1  j-i  A 

f.y 

1  j-1,k 

l  J  J 

l  j  J 

+  | 

I  +  1 

[g  V 

1 

l  J  J 

jA+l  - 

i  J  J 

1  J,k+1 

fris) 

jA-l  1 

[g  V] 

|  ] 

l  J  J 

i  J  J 

Nov  replacing  the  metric  derivatives  with  their 
finite  difference  expressions,  equations  (7) 
through  (11),  and  writing  out  only  the  continuity 
equation  for  clarity,  we  have 


f  iP/J.) n  1  J 

(  In  | 

Hi*i,kiyj-i^i  -  yj+i.k-ij 

■  H"+i.kh+i'k+i  ■  x2+,'k-’) 

+  [‘  H",k+i(yi+1'k+1  -  vj-i.k.i) 

+  Hj.k+il*}+1'k+1  * 


-  -  (eu)j,k-i(yj+i.v-i  -  yj-i,k-i) 

+  HU-b*'*-'  ■  x3-i'k-i)]J 


-  -  (RC)jfk  .  (15 

This  form  of  the  continuity  equation  is  identical 
to  that  which  would  be  obtained  by  an  integration 
of  the  control  volume  fluid  equations  around  the 
cell  area,  assuming  no  variation  of  (pu)  or  (pv) 
across  a  cell  face. 

Again  using  only  the  continuity  equation  for 
clarity, 

^  ^pdxdy  =  ix  +  V  **  (16 

area  control  surface 


where  i„,  i„  are  unit  vectors  in  original 

x  y 

Cartesian  coordinates, 

n  is  an  outward  normal  vector  to 

the  cell  face,  and 

dl  is  length  along  the  face. 

Equation  (16),  when  expanded,  is  identical  to 
equation  (14),  and  shows  that  the  proper 
interpretation  of  a  metric  derivative  is  that  the 
ratio  of  that  derivative  to  the  Jacobian  is  the 
projected  area  of  the  face.  Thus  the  solution  of  a 
finite  difference  problem  in  strong  conservation 
law  form  in  a  transformed  space  is  equivalent  to 
the  solution  of  a  flux  summation  algorithm 
expressed  in  physical  space. 

The  flux  sum  of  equation  (13)  could  be 
improved  by  selecting  a  more  accurate 
approximation  for  the  integral  of  (pu)  or  (pv) 
over  a  face  than  simply  its  value  at  the  center 
of  the  face.  For  example,  a  piecewise  continuous 
interpretation  of  (pu)  is  equivalent  to 

J*(eu  ix)’K  dl  =  (vk(Pu>j,k]  4k  Vj,k 

+  |vk<ou)j(Vj  «k  yj(k  .  (17) 

An  important  property  of  equation  (14)  is  that 
a  uniform  solution  [(pu)j^k  ,  (pv)j^k  constant 
over  space]  remains  identically  uniform.  This 
property  is  generally  referred  to  as  area  conserva¬ 
tion  in  finite  element  analysis  and  its  importance 
to  maintaining  constant  freestream  solutions  was 
first  pointed  out  by  Steger  (7).  A  recent  paper  by 
Hindman  [9]  has  also  emphasized  this  property. 
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COMPARI  SON  OF  FLUX  BALANCE  FORM  OF  SEVERAL 

TIME -MARCHING  ALGORITHMS  FOR  INVISCID  FLOW  The  flux  balance  operator  is: 


Centered  Difference  Algorithms 


METHOD  A.  Implicit  Time-Marching  Algorithm  of 
Beam  and  Warming  [1]: 

The  first  method  to  be  considered  is  the 
time  marching  algorithm  introduced  by  Beam  and 
Warming  which  is  unconditional ly  stable,  in  the 
linear  problem,  for  two  dimensional  flows: 

(i  +  At  (i  ♦  At  L*)(G£J  -  «j,k)  ' 

“  -  At^«°  rc-j.k  +  4k  ^j.k  "  lBWAV  G  j  ,kj 
=  -  At{RCj^  -  Lg^y  uj,k^  '  *  11 

where  L^  and  L^  are  finite  difference 
operators  whose  forvr*  determines  the  order  of 
accuracy  of  the  time-marching  algorithm  and  L gvJAV 
is  a  finite  difference  smoothing  operator. 


Solution  Procedure 


Step  1 : 

(l  ♦  L{j  *  "  At(RC:,k  “  LBWAVUj,k]  (19) 

Step  2: 

(i  +  st  Ln)  4Gj,k  -  aG*, 


i*  * 


(20) 


Step  3: 


~n+ 1  'sn  ''n 

uj,k  -  uj,k  +  Auj,k 


(21) 


METHOD  B.  Modified  Beam  and  Warming  Algorithm 

In  order  to  improve  the  accuracy  of  the  flux 
balance  operator  in  skewed  mesh  systems,  a  variety 
of  modified  Beam  and  Warming  algorithms  have  been 
developed.  The  simplest  modified  operator,  which 
corresponds  to  several  other  algorithms,  is  to 
balance  the  flux  on  the  cell  showv. -below. 


j  ,k+1 

j+1 ,  k+ 1 

i-i 

M 

k+* 

i 

y 

I  2 

X 

j-1,k 

* 

J-i' 

k-A 

i 

is  ”  i 
_ i 

,.1  V_i 
i'K  1 

j+lrk 

jrX-1 

j-M,k-1 

Fioure  4.  Half-spacing  resh  cell  in  computational 
space. 


RCHj,k  ■  4j/2<FCH)"(k  +  i°/2(GCH)j>k  ,  (22) 

where 

(FCH)"tJA  -  (kj  Fj,k)(“j  «k/2  yj,k) 

-  (vj  Gj,k][«j  4k/2  x  j  ,kj 

(GCH)"(k±i  7  -  (“k  Fj#k)  (wk  {j/2  y j ,kj 

+  (“k  Gj,k]  (uk  4j/2  xj,k]  < 

and  the  algorithm  becomes 

(I  -  Atj(kL^)(I  +  Atj>kL^)(u"*k  -  G"(k] 

=  ”  4tj »k ( (RCH> j ,k  ”  ^BWAV  uj,k]  •  (23' 


METHOD  C.  Explicit  Time  Integration  Scheme  of 
Jameson  (3]  and  Rizzi  [4]  : 

The  Runae-Kutta  time  stepping  approach  can  be 
expressed  as: 


dU 


dt 


-  -  RCH5,k  -  LrxAV  G"(k  -  -  pG",k 


Step  1 : 


u<0)  -  «j,k 


Step  2: 


G(i)  =  G'°> 


{  2Jj. 

G(2)  .  G<0)-  f— )  pG(,) 

l  2Jj,x 


Step  3: 


Step  4 : 


M3)  MO)  .  .  p  ( 2 ) 
U  =  U  -  Atj^  PU 


Step  5: 


G(4)  .  G(0> 


(24) 


(25) 


(26) 


(27) 


( 2B ) 


-  M  (pG(0).  2pGM>.  2pG(2)+  pG<3,|  (29) 
V  6Jj,kV 


Step  6 

Gn+1  -  G{4) 

Uj,k  “  u 


(30) 
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This  method  is  usually  run  with  the 
smoothing  operator  in  each  step  formed  as 

kjuCAV  U  rather  than  ^RKAV  ^  * 

If  steps  2,  3,  4,-5  each  converge  independent¬ 
ly,  then  the  steady  state  solution  is  identical  to 
that  of  the  modified  Beam  and  Warming,  method  B,  if 
**BWAV  *  LRKAV'  an<3  solution  is  independent 

Of 

METHOD  D.  Flux  Balance  Form  of  MacConnack's 
Method  f 2] : 

A  flux  balance  form  of  MacCormack* s  method 
may  be  expressed  as: 


Step  1,  Predictor: 


Step  3  may  be  rewritten  as 
Ujtk  “  Uj,k 


Us 


(<Fj+l,k  +  Fj,k><»j  4k/2  Yj,k> 

*  tGj+1,k  +  Gj,k}<yj  4  k/2  xj,k>] 

“  (tFj,k  +  Fj-1,kHuj  *k/2  yj#k) 

"  <Gj,k  +  G3-1,kHyj  4k/2  xj,k>] 

+  (“  {Fj,k+1  +  Fj,k}<yk  4j/2  Vj,k> 


(  n  n  \ 

*  °j,k  "  Atj,k 

[Fj+1,kPLy1  *  Gj+1,kPLx1 

+  (Gj,k+1  +  1 

-  (Fj,kPLy3  '  G 

j,kPLx3  j 

-  (-  <Fj,k  + 

+  ("  Fj,k+1PLy4 

+  Gj,k+1PLx4| 

+  «Gj,k  +  Gj 

n  'S 

•  [-  Fj,kPLy2  + 

Gj,kPLx2  1 

“  ^MAV  ^Uj,k 

(34) 


“  LmAV  uj ,kj 
where  the  projected  lengths  are 
PLy1  -  Uj  «k/ 2  Yj,k  ■  Vj+i,k+i  •  yj+l.k-} 


(31)  This  form  of  MacCormack's  method  is  area 

preserving  for  each  step  (31  or  32)  and  offers 
the  possibility  of  steady  state  solutions  which 
are  independent  of  At.  When  the  predictor  and 
corrector  step  both  converge,  which  we  generally 
find  to  be  the  case,  then  the  steady  state  form 
of  equation  (34)  is: 


PLX4  k  4j/2  xj,k  “  xj*i,k+^  * 


RCH 


j ,  k  “  h^AV  uj,k  “  0  • 


(35) 


as  with  the  modified  Beam  and  Warming  and  fourth 
order  Runge-Kutta  schemes,  methods  B  and  C. 


Step  i 

l,  Corrector : 

^  #  T/ 

•  < 

uj,k  -  Uj.k  -  itj.kN 

Fj,kPLy1  “  G 

- 1 

[F}-1,kPLy3  " 

Gj-1,kPLx3  ) 

♦  i 

(-  Fj>kPI"y4  + 

Gj ,kPLx4j 

- 1 

(*  Fj,k-1pS-2 
-.  1 

+  Gj,k-1PLx2 

■  lmav  uj ,k 

* 

Step  3: 

Sft  ■  7  (Sj-lc  +  "i\) 


(32) 


(33) 


Comments  on  Centered  Difference  Aloorithms 


For  each  of  the  methods  considered.  A,  B,  C  and 
D,  the  steady  state  solution  desired  is  that  the 
centered  differenced  flux  balancing  residual 
(RC)^y,  or  (  RCH )  -j  is  zero.  For  all 
these  schemes  at  interior  points,  this  criterion  is 
possible  only  when  the  artificial  smoothing 
operators  are  neglected.  The  order  of  accuracy  of 
all  these  methods  is  the  same  and  each  should  be 
expected  to  produce  nearly  the  same  steady  state 
solution  except  for  differences  in  smoothing 
requirements  or  boundary  conditions.  As  we  shall 
see  in  the  test  examples,  this  conclusion  can  be 
demonstrated  by  numerical  experiments. 


Nor.-Centered  Difference  Algorithms 


M FT H OD JE .  MacCormack 1 s  Method  on  a  Non-Centered 
Cell  Basis: 

A  version  of  MacCormack • s  method  which  Joes  a 
flux  sum  on  the  cell  structure  shown  in  Figure  5 
has  been  developed  by  Tong  (10).  This  form  is  a 
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where 


two  or  three  dimensional  extension  of  the  strong 
conservation  law  version  of  MacCormack's  method 
proposed  by  Hindman  [9] .  This  scheme  uses  the  same 
basic  nine  point  cell  as  in  Figure  3  or  4,  but  it 
computes  a  flux  sum  over  each  of  the  four  interior 
cells. 

Step  1 : 

uj,k  “  uj,k  ~  j  ,k  RMC )  ,k+j  “  lMCAV  uj,k 


Step  2: 


uj,k  =  Uj,k  “  At j / k|^ RMC J j- j #k—J  "  lMCAV  uj,kj 
Step  3; 

-  j  («5.k  ♦  «3Tk) 


uj!k  ■  uj,k 


Ati  k 

~2  |<RMC)*j+i,k+i 


+  (RMC)  j--i#k-i  *  LMCAV^uj,k  4  uj,k* 
where 

( RMC)n+ 1  fc+1  -  fl  l”tiS.k 

3**'ki*  Aji*,k±i 

(FM_)j,k  =  F  j , k^ kV j ,  k  ~  Gj,k*kxj,k  ' 


~  “  Fj#k6jyj,k 


*3  jxj#>< 


(36) 


(37) 


(38) 


(39) 


(40) 


(41  ) 


(42) 


*j±1/2,ki1/2  t^,e  area  4-point  cell. 

This  form  of  MacCormack's  method  also  obeys 
geometric  area  conservation  on  the  4-point  cells 
with  the  possibility  of  a  steady  state  solution 
that  is  independent  of  At.  Both  the  predictor 
and  corrector  steps  can  individually  converge. 
The  predictor  and  corrector  steps  express  flux 
conservation  on  the  4-point  cells  and,  if  each 
( RMC )  j  ,  k  converges  to  zero,  then  we  have  global 
flux  conservation. 

METHOD  F.  Ni 1 s  Lax-Wend ro f f  Method  [5]  : 

Ni’s  method  may  be  expressed  in  a  one  step 
form  as: 


nn+1 

uj,k 


*  uj , k  *  Atj,k<Hj/2  y k/2  (RNI)j,k} 
+  (u°/2«j/2[ 


f  bmi)  ( p 

Jill 

{  J  V 

< 

n 


f  3GNl)  ( RNl) 

n  1 

\  SU  J  l  A  )\ 

1  j.kj 

4  ^NIAV  uj,k 


(43) 


(rni) 


(44) 


(45) 


^j±i,k±i 

1  (“k  Fj,k)(4k  yj,k]  -  (“k  Gj,k)(6k  *j(k] 

t®”)***.* 

=  •  (“j  rj/k)  [4j  Vj.k)  +  [A  Gj,k)({j  xj,k)  <46) 

This  representation  of  Ni  scheme,  which  is 
not  the  computational  form,  is  a  Lax-Wendroff 
type  method  which  does  a  flux  sum  on  4-point 
cells.  A  possible  solution  to  the  steady  state 

form  of  equation  (43)  is  that  each  (  RNI )  ^ /2  ,  )<±  1/2 

becomes  zero  in  the  absence  of  artificial 
smoothing . 

Comments  on  Non-Centered  Schemes 

Both  the  Tong  and  the  Ni  non-centered  schemes 
compute  flux  balances  on  the  4-point  cells  and 
each  solution  can  be  consistent  with  these 

balances,  ( RMC) j+ i/2 ,k± 1/2  or  ^ RNI ^ j± 1/2 , ki 1/2  ' 

becoming  zero  in  the  absence  of  smoothing.  Thus  the 
Tong  and  Ni  schemes  should  produce  approximately 
the  same  steady-state  solution,  if  the  same 
boundary  conditions  are  applied.  The  two  steady 
state  flux  balance  operators  are  not  identical.  We 
would  expect  Ni’s  scheme  to  produce  a  more  accurate 
result  in  the  absence  of  artificial  smoothing,  but 
his  method  may  require  more  smoothing  than  Tong's 
method . 

INVISCID  TEST  EXAMPLES 

Three  configurations  were  chosen  as  test 
examples.  The  first  case  was  the  bump  in  a  channel 
problem  introduced  by  Ni  (5);  the  second  was  a 
supersonic  nozzle  with  a  rapid  isentropic 
expansion;  and  the  third  was  a  supercritical  stator 
designed  by  Sanz  (11).  For  each  problem,  an  exact 
steady  state  solution  would  be  isentropic  with 
constant  stagnation  pressure,  stagnation 
temperature  and  mass  flow  rate.  Methods  are 
compared  on  the  basis  of  predicted  Mach  number 
distributions  and  stagnation  pressure  loss  over  the 
domain.  For  all  methods  on  all  problems,  the  mass 
flow  rate  was  conserved  to  within  0.3%. 

Boundary  Conditions 

Inflow/outflow  boundary  conditions  were  of  the 
non-reflecting  characteristic  or  extrapolation 
type,  but  implementation  details  differed  greatly 
from  method  to  method.  We  believe  that  the  results 
presented  are  as  free  as  possible  from  contamina¬ 
tion  due  to  these  conditions. 
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It  was  impossible  to  maintain  exactly  the  same 
hard  wall  boundary  conditions  between  the  schemes. 
Methods  A  and  B#  which  were  developed  strictly  for 
viscous  flow  by  the  authors#  used  only  a  simple 
pressure  and  parallel  velocity  extrapolation  for 
the  inviscid  case.  These  conditions  were  sufficient 
for  the  first  two  cases  but  not  the  supercritical 
stator  case.  The  proper  conditions  for  Ni's  method 
are  not  clear.  The  first  two  examples  for  Ni's 
method  predicted  surface  quantities  from  the 
interior  point  scheme  and  then  corrected  these 
quantities  by  a  simple  one-dimensional  wave  to 
maintain  flow  tangency. 

Methods  C,  D#  and  E  used  boundary  procedures 
which  consistently  satisfied 

normal  momentum 

1 

lE  -  (47) 

3n  R 

ill 

where  q  *  u  +  v 

and  R  is  the  radius  of  curvature  of  the  surface 
streamwise  momentum 

ip  *  -  0  3 ( q  /2 ).  (48) 

3s  3s 

where  s  is  the  arc  length  along  surface 
flow  tangency 

v/u=  6jy  /  6  jX  (49) 

at  each  point  on  the  boundary.  To  satisify  all 
three  conditions  at  each  point  on  the  boundary 
requires  an  iterative  evaluation  of  the  boundary 
pressures  at  each  time-marching  step.  This 
procedure  is  costly  but  is  essential  for 
simulations  with  minimum  stagnation  pressure  error. 


Art i f i ci al  Smoothing  Operators 


For  methods  A#  B#  D,  and  F,  the  smoothing 
operators  were  generally  as  proposed  by  their 
authors  with  the  smoothing  coefficients  "tuned-up” 
for  each  calculation.  For  methods  C,  D,  and  E,  the 
calculation  procedure  was: 


[I 

-  RSFB] Un 

(50) 

tl 

-  SX)U* 

(51) 

(I 

♦  Sy)U** 

(52) 

where  RSFB  is  the  flux  balance  operator  and 

S„  #  Sv  are  the  smoothing  operators 
x  y 

Sxu  -  f[  Uj+1  +  U i - 2U j] 
f  [l  _fJ_‘ 

1  +  f2  [_  Ax/Ax^in 

f  1  -  0. 1  f2  «  200. 


(53) 

(54) 


For  methods  C  and  D#  these  smoothing  operators  were 
applied  after  each  step#  but  for  method  E  they  were 
applied  only  after  the  last  step. 

The  steady  state  solutions  are  not  independent 
of  the  smoothing  operators#  but  every  effort  was 
made  to  reduce  their  effect. 


CONVERGENCE  CRITERION 

For  all  of  the  methods  except  Tong's#  the 
convergence  criterion  was  easy  to  establish.  The 
computed  solutions  were  run  until  the  steady  state 
flux  balance  operator  plus  the  smoothing  operator 
converged  to  a  specified  value#  1.0  X  10"^ 
maximum  at  any  point  and  1.0  x  10"^  root  mean 
square  value. 

In  Tong's  method#  as  with  any  MacCormack  type 
method#  the  order  of  differencing  in  predictor  and 
corrector  steps  is  not  unique.  The  double  sawtooth 
solutions  common  to  central  difference  schemes  can 
be  reduced  by  cyclically  permuting  the  operator 
sequence.  The  disadvantage  of  this  process  plus  the 
smoothing  sequence  adopted  is  that  we  lose  the 
ability  to  drive  a  steady  state  residual  below 
about  1.0  x  10~4  maximum  local  value  and  about 
1.0  x  10~5  root  mean  square  value.  It  is  thus 
difficult  to  determine  when  to  stop  iterating  on  a 
solution#  and  we  must  supply  some  heuristic 
principle  to  determine  when  to  stop.  The  advantage 
of  this  process  is  that  the  stagnation  pressure 
errors  in  the  "converged”  solution  are  greatly 
reduced. 


Ni  Bump  Test  Case,  Subsonic  Flow 

All  of  the  methods  produced  acceptably  accurate 
solutions  for  the  channel  bump  problem  especially 
in  terms  of  symmetry  and  Mach  number  predictions. 
The  grid  used  contained  65X17  points  and  is 
illustrated  in  Figure  5  along  with  a  typical  Mach 
number  contour  plot.  The  predicted  top  and  bottom 
wall  Mach  numbers  for  each  method  are  shown  in 
Figure  6.  The  stagnation  pressure  loss  for  each  of 
the  methods  is  shown  in  Figure  7.  All  of  the 
centered  schemes  produced  essentially  the  same 
pressure  loss  of  about  2%  max  APt/Pt.  The 
non-centered  schemes  produced  essentially  no  loss. 
On  balance#  all  the  schemes  performed  adequately  on 
the  bump  case. 


Supersonic  Nozzle  Test  Case 

The  supersonic  nozzle  test  case  geometry 
and  grid  is  shown  in  Figure  8.  This  grid  has  65X33 
points.  The  inflow  Mach  number  is  2.0#  and  an 
idealized  wave  diagram  is  also  shown  in  Figure  8. 
The  wave  diagram  shows  that  the  wall  outflow  Mach 
number  should  be  2.83#*  and  that  the  Mach  number 
should  be  uniform  in  region  1.  The  idealized  wave 
diagram  does  not  account  for  curvature  of  the 
characteristics  as  the  wave  system  from  the  lower 
wall  interacts  with  the  upper  wall  wave  system. 
Under  this  assumption  the  exit  centerline  Mach 
number  should  be  2.625.  The  results  for  each  of 
the  schemes  is  shown  in  Figure  9  in  terms  of  Mach 
number  contours.  The  correct  wave  system  is 
adequately  predicted  by  all  schemes.  The 
stagnation  pressure  loss  contours  for  this  test 
case  are  all  shown  in  Figure  10.  The  centered 
schemes  all  have  unacceptable  losses  in  the  range 
of  3%  to  6%.  Ni's  scheme  has  a  small  loss  while 
the  non-centered  MacCormack  has  virtually  no  loss. 
It  is  likely  that  a  better  hard  wall  boundary 
condition  treatment  with  Ni's  method  would 
eliminate  the  loss  shown. 
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The  agreement  between  the  4  centered 
schemes  is  quite  remarkable;  in  fact,  as  would  be 
expected  from  the  analysis,  all  schemes  produce 
nearly  the  same  steady  state  solutions.  The 
nor.-centered  MacCormack  and  Ni's  scheme  produce 
quite  similar  results. 

Supercritical  Stator  Test  Case 

A  much  more  demanding  geometry  is  the 
supercritical  stator  design  presented  by  Sant  111). 
The  grid  for  this  test  case  is  shown  in  Figure  11. 

A  simple  sheared  grid  was  used  for  simplicity  which 
required  a  fine  grid  near  the  blade  leading  edge 
for  accurate  simulations. 

The  surface  Mach  number  predictions  for 
the  supercritical  stator  using  the  Tone’s  non- 
centered  MacCormack  type  method  are  shown  in  Figure 
12.  The  agreement  is  satisfactory  except  at  8% 
chord  on  the  suction  surface.  Slight  variations  in 
the  specified  downstream  pressure  specified 
produced  overshoots  or  undershoots  in  this  region. 
The  sample  calculation  is  for  the  downstream 
pressure  specified  by  Sanz.  We  are  simply  unsure 
of  the  reason  for  this  discrepancy.  The  stagnation 
pressure  loss  for  this  calculation  is  shown  in 
Figure  ’3.  An  unacceptably  high  loss  of  14%  exists 
locally  at  the  stagnation  point  and  a  maximum  loss 
of  4%  exists  on  the  suction  surface. 

The  centered  MacCormack  scheme  performed  nearly 
the  same  on  the  Sanz  blade  as  is  shown  in  Figure 
14.  The  stagnation  pressure  loss  at  the  leading 
edge  is  locally  14%  and  a  maximum  of  4%  in  the 
remainder  of  the  domain. 


VISCOUS  TURBINE  CALCULATION  EXAMPLE 

While  the  results  from  the  supercritical 
stator  geometry  illustrate  that  accurate  Euler 
equation  solutions  can  be  obtained,  the  situation 
is  less  clear  for  viscous  flow  problems.  To 
illustrate  our  level  of  computational  abilities  for 
these  flows,  a  calculation  of  the  flow  in  a  high 
speed,  high  turning  turbine  cascade  is  presented. 
The  numerical  method  used  is  similar  to  the 
approximate  factorization  method  of  Beam-Warming 
[1]  except  that  the  flux  balance  operator  uses 
trapezoidal  rule  integration,  as  in  equation  (17), 
on  the  computational  cell  of  Figure  4.  Local 
values  of  At  are  used  to  produce  a  constant  local 
CFL  number  of  about  4.  Details  of  this  alaorithm 
will  be  available  in  reference  [12].  Solutions  for 
turbine  geometries  generally  require  about  500 
iterations  to  converge  to  machine  accuracy,  6 
decimal  digits. 

The  experimental  results  and  geometric  details 
for  the  cascade  were  made  available  in  reference 
[13],  and  the  geometry  is  shown  in  Figure  15.  The 
turning  is  126  degrees  with  an  outflow  Mach 
number  of  0.75.  The  design  Reynolds  number  is 
500000.  All  calculations  assumed  laminar  flow. 

The  computational  results  are  presented  in 
Figures  16  and  17.  The  first  figure  compares 
predicted  blade  surface  pressure  to  the 
experimental  values  for  a  laminar  calculation  at 
design  Reynolds  number.  The  grid  chosen  has  100x29 
points  and  is  nearly  orthogonal  at  eacn  node  point. 


Predicted  results  are  quite  sensitive  to  the  grid 
resolution.  Reasonably  accurate  surface  pressures 
were  predicted  everywhere  including  the  stagnation 
point  and  the  trailing  edge.  The  maximum  stagnation 
pressure  error  outside  the  boundary  layer  was  about 
1%.  We  are  as  yet  unable  to  explain  the  discrepancy 
at  70%  chord  on  the  suction  surface  although  we 
believe  it  to  be  linked  to  the  trailing  edge  flow 
field. 

The  predicted  trailing  edge  flow  field  is 
illustrated  in  Figure  17  in  terms  of  velocity 
direction  vectors.  This  calculation  shows  two  small 
separation  zones  comparable  in  size  to  the  trailing 
edge  radius.  While  the  predicted  flow  field  at  the 
trailing  edge  seems  reasonable,  we  have  no  way  to 
verify  its  accuracy. 


CONCLUSIONS  AND  DISCUSSION 


The  presentation  of 
schemes  in  terms  of  a  flux 
and  the  inviscid  computati 
the  six  schemes  analyzed 
steady  state  solution  prr 
computational  details  of 
different,  the  steady  sta 
central  difference  scheme 
Computational  details  of 
are  also  quite  different 
similar  steady  state  sol. 


the  time-marching 
balancing  interpretation 
•al  examples  shows  that 
■nuite  similar  in  their 
ties.  While  the 
schemes  are  very 
’utions  of  the 
-e  nearly  identical, 
ion-centered  schemes 
ev  again  produce 


A  consistent  result  throughout  the 
evaluations  was  that  the  non-centered  schemes 
produced  lower  stagnation  pressure  errors  than  did 
the  centered  schemes.  Since  entropy  conservation  or 
stagnation  pressure  conservation  is  not  a 
conservation  property  of  the  finite  difference 
equations,  we  expect  these  errors  to  be  a  function 
of  the  truncation  errors  present  in  the  solution. 
The  truncation  error  appears  as  an  entropy 
generation,  and  the  stagnation  pressure  error  is 
the  entropy  error  raised  to  the  t  /  ( 3r  3 )  power. 
Stagnation  pressure  loss  is  a  very  sensitive 
measure  of  the  entropy  generation  in  a  scheme. 


The  grid  sizes  for  the  inviscid  problems  were 
chosen  to  produce  small  stagnation  pressure  errors 
for  the  non-centered  schemes,  and  the  centered 
schemes  nearly  always  produced  larger  errors  on  the 
same  grid.  We  attribute  this  result  either  to  the 
space  averaging  properties  of  the  solution 
operators  or  to  artificial  smoothing  to  eliminate 
double  sawtooth  errors. 


Accurate  Fuler  flow  simulations  are  possible 
with  the  non-centered  schemes  as  was  shown  by  the 
supe rcri t i cal  stator  example.  Both  methods 
reproduced  the  shock-free  solution  with  a  minimum 
of  effort,  and  both  produced  only  a  small 
stagnation  pressure  loss  on  the  test  grid.  With  a 
locally  refined,  orthogonal  grid,  the  total  number 
of  grid  points  required  can  be  reduced  to  a  nore 
manageable  number.  Our  ability  to  compute  accurate 
viscous  flow  solutions  is  increasing  rapidly.  The 
turbine  cascade  solution  is  for  a  moderately 
difficult  geometry,  and  the  computational  results 
are  reasonable.  Much  more  work  must  be  done  in  this 
area  before  we  can  compute  viscous  solutions  as 
accurately  as  for  inviscid  situations. 
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Ficure  6a.  Surface  .Mach  number,  Ni  bump, 
Beam  Warning  method 


[13]  private  communication  from  Eolls-Royce 

Figure  6b.  Surface  Mach  number,  Ni  tunp, 
modified  Be an-Warrinc  method 
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Figure  16  Blade  surface  static  pressure 
comparison  for  turbine  stator 


Figure  17  Velocity  vector  direction  plot  for 
turbine  stator  trailing  edge  region 


